module SnowHydrologyMod

!-----------------------------------------------------------------------
!BOP
!
! !MODULE: SnowHydrologyMod
!
! !DESCRIPTION:
! Calculate snow hydrology.
!
! !USES:
  use shr_kind_mod, only: r8 => shr_kind_r8
  use clm_varpar  , only : nlevsno
!
! !PUBLIC TYPES:
  implicit none
  save
!
! !PUBLIC MEMBER FUNCTIONS:
  public :: SnowWater         ! Change of snow mass and the snow water onto soil
  public :: SnowCompaction    ! Change in snow layer thickness due to compaction
  public :: CombineSnowLayers ! Combine snow layers less than a min thickness
  public :: DivideSnowLayers  ! Subdivide snow layers if they exceed maximum thickness
  public :: BuildSnowFilter   ! Construct snow/no-snow filters
!
! !PRIVATE MEMBER FUNCTIONS:
  private :: Combo            ! Returns the combined variables: dz, t, wliq, wice.
!
! !REVISION HISTORY:
! Created by Mariana Vertenstein
!
!EOP
!-----------------------------------------------------------------------

contains

!-----------------------------------------------------------------------
!BOP
!
! !IROUTINE: SnowWater
!
! !INTERFACE:
  subroutine SnowWater(lbc, ubc, num_snowc, filter_snowc, &
                       num_nosnowc, filter_nosnowc)
!
! !DESCRIPTION:
! Evaluate the change of snow mass and the snow water onto soil.
! Water flow within snow is computed by an explicit and non-physical
! based scheme, which permits a part of liquid water over the holding
! capacity (a tentative value is used, i.e. equal to 0.033*porosity) to
! percolate into the underlying layer.  Except for cases where the
! porosity of one of the two neighboring layers is less than 0.05, zero
! flow is assumed. The water flow out of the bottom of the snow pack will
! participate as the input of the soil water and runoff.  This subroutine
! uses a filter for columns containing snow which must be constructed prior
! to being called.
!
! !USES:
    use clmtype
    use clm_varcon  , only : denh2o, denice, wimp, ssi
    use clm_time_manager, only : get_step_size
    use clm_atmlnd        , only : clm_a2l
    use SNICARMod         , only : scvng_fct_mlt_bcphi, scvng_fct_mlt_bcpho, &
                                   scvng_fct_mlt_ocphi, scvng_fct_mlt_ocpho, &
                                   scvng_fct_mlt_dst1,  scvng_fct_mlt_dst2,  &
                                   scvng_fct_mlt_dst3,  scvng_fct_mlt_dst4
!
! !ARGUMENTS:
    implicit none
    integer, intent(in) :: lbc, ubc                    ! column bounds
    integer, intent(in) :: num_snowc                   ! number of snow points in column filter
    integer, intent(in) :: filter_snowc(ubc-lbc+1)     ! column filter for snow points
    integer, intent(in) :: num_nosnowc                 ! number of non-snow points in column filter
    integer, intent(in) :: filter_nosnowc(ubc-lbc+1)   ! column filter for non-snow points
!
! !CALLED FROM:
!
! !REVISION HISTORY:
! 15 September 1999: Yongjiu Dai; Initial code
! 15 December 1999:  Paul Houser and Jon Radakovich; F90 Revision
! 15 November 2000: Mariana Vertenstein
! 2/26/02, Peter Thornton: Migrated to new data structures.
! 03/28/08, Mark Flanner: Added aerosol deposition and flushing with meltwater
!
! !LOCAL VARIABLES:
!
! local pointers to implicit in arguments
!
    integer , pointer :: snl(:)              !number of snow layers
    logical , pointer :: do_capsnow(:)       !true => do snow capping
    real(r8), pointer :: qflx_snomelt(:)     !snow melt (mm H2O /s)
    real(r8), pointer :: qflx_rain_grnd(:)   !rain on ground after interception (mm H2O/s) [+]
    real(r8), pointer :: qflx_sub_snow(:)    !sublimation rate from snow pack (mm H2O /s) [+]
    real(r8), pointer :: qflx_evap_grnd(:)   !ground surface evaporation rate (mm H2O/s) [+]
    real(r8), pointer :: qflx_dew_snow(:)    !surface dew added to snow pack (mm H2O /s) [+]
    real(r8), pointer :: qflx_dew_grnd(:)    !ground surface dew formation (mm H2O /s) [+]
    real(r8), pointer :: dz(:,:)             !layer depth (m)
!
! local pointers to implicit out arguments
!
    real(r8), pointer :: qflx_top_soil(:)     !net water input into soil from top (mm/s)
!
! local pointers to implicit inout arguments
!
    real(r8), pointer :: h2osoi_ice(:,:)     !ice lens (kg/m2)
    real(r8), pointer :: h2osoi_liq(:,:)     !liquid water (kg/m2)
    integer , pointer :: cgridcell(:)        ! columns's gridcell (col) 
    real(r8), pointer :: mss_bcphi(:,:)      ! hydrophillic BC mass in snow (col,lyr) [kg]
    real(r8), pointer :: mss_bcpho(:,:)      ! hydrophobic BC mass in snow (col,lyr) [kg]
    real(r8), pointer :: mss_ocphi(:,:)      ! hydrophillic OC mass in snow (col,lyr) [kg]
    real(r8), pointer :: mss_ocpho(:,:)      ! hydrophobic OC mass in snow (col,lyr) [kg]
    real(r8), pointer :: mss_dst1(:,:)       ! mass of dust species 1 in snow (col,lyr) [kg]
    real(r8), pointer :: mss_dst2(:,:)       ! mass of dust species 2 in snow (col,lyr) [kg]
    real(r8), pointer :: mss_dst3(:,:)       ! mass of dust species 3 in snow (col,lyr) [kg]
    real(r8), pointer :: mss_dst4(:,:)       ! mass of dust species 4 in snow (col,lyr) [kg]
    real(r8), pointer :: flx_bc_dep_dry(:)   ! dry BC deposition (col) [kg m-2 s-1]
    real(r8), pointer :: flx_bc_dep_wet(:)   ! wet BC deposition (col) [kg m-2 s-1]
    real(r8), pointer :: flx_bc_dep(:)       ! total BC deposition (col) [kg m-2 s-1]
    real(r8), pointer :: flx_bc_dep_pho(:)   ! hydrophobic BC deposition (col) [kg m-1 s-1]
    real(r8), pointer :: flx_bc_dep_phi(:)   ! hydrophillic BC deposition (col) [kg m-1 s-1]
    real(r8), pointer :: flx_oc_dep_dry(:)   ! dry OC deposition (col) [kg m-2 s-1]
    real(r8), pointer :: flx_oc_dep_wet(:)   ! wet OC deposition (col) [kg m-2 s-1]
    real(r8), pointer :: flx_oc_dep(:)       ! total OC deposition (col) [kg m-2 s-1]
    real(r8), pointer :: flx_oc_dep_pho(:)   ! hydrophobic OC deposition (col) [kg m-1 s-1]
    real(r8), pointer :: flx_oc_dep_phi(:)   ! hydrophillic OC deposition (col) [kg m-1 s-1]
    real(r8), pointer :: flx_dst_dep_dry1(:) ! dry dust (species 1) deposition (col) [kg m-2 s-1]
    real(r8), pointer :: flx_dst_dep_wet1(:) ! wet dust (species 1) deposition (col) [kg m-2 s-1]
    real(r8), pointer :: flx_dst_dep_dry2(:) ! dry dust (species 2) deposition (col) [kg m-2 s-1]
    real(r8), pointer :: flx_dst_dep_wet2(:) ! wet dust (species 2) deposition (col) [kg m-2 s-1]
    real(r8), pointer :: flx_dst_dep_dry3(:) ! dry dust (species 3) deposition (col) [kg m-2 s-1]
    real(r8), pointer :: flx_dst_dep_wet3(:) ! wet dust (species 3) deposition (col) [kg m-2 s-1]
    real(r8), pointer :: flx_dst_dep_dry4(:) ! dry dust (species 4) deposition (col) [kg m-2 s-1]
    real(r8), pointer :: flx_dst_dep_wet4(:) ! wet dust (species 4) deposition (col) [kg m-2 s-1]
    real(r8), pointer :: flx_dst_dep(:)      ! total dust deposition (col) [kg m-2 s-1]
    real(r8), pointer :: forc_aer(:,:)       ! aerosol deposition from atmosphere model (grd,aer) [kg m-1 s-1]
!
!
! !OTHER LOCAL VARIABLES:
!EOP
!
    integer  :: c, j, fc                           !do loop/array indices
    real(r8) :: dtime                              !land model time step (sec)
    real(r8) :: qin(lbc:ubc)                       !water flow into the elmement (mm/s) 
    real(r8) :: qout(lbc:ubc)                      !water flow out of the elmement (mm/s)
    real(r8) :: wgdif                              !ice mass after minus sublimation
    real(r8) :: vol_liq(lbc:ubc,-nlevsno+1:0)      !partial volume of liquid water in layer
    real(r8) :: vol_ice(lbc:ubc,-nlevsno+1:0)      !partial volume of ice lens in layer
    real(r8) :: eff_porosity(lbc:ubc,-nlevsno+1:0) !effective porosity = porosity - vol_ice
    integer  :: g                                 ! gridcell loop index
    real(r8) :: qin_bc_phi(lbc:ubc)               ! flux of hydrophilic BC into layer [kg]
    real(r8) :: qout_bc_phi(lbc:ubc)              ! flux of hydrophilic BC out of layer [kg]
    real(r8) :: qin_bc_pho(lbc:ubc)               ! flux of hydrophobic BC into layer [kg]
    real(r8) :: qout_bc_pho(lbc:ubc)              ! flux of hydrophobic BC out of layer [kg]
    real(r8) :: qin_oc_phi(lbc:ubc)               ! flux of hydrophilic OC into layer [kg]
    real(r8) :: qout_oc_phi(lbc:ubc)              ! flux of hydrophilic OC out of layer [kg]
    real(r8) :: qin_oc_pho(lbc:ubc)               ! flux of hydrophobic OC into layer [kg]
    real(r8) :: qout_oc_pho(lbc:ubc)              ! flux of hydrophobic OC out of layer [kg]
    real(r8) :: qin_dst1(lbc:ubc)                 ! flux of dust species 1 into layer [kg]
    real(r8) :: qout_dst1(lbc:ubc)                ! flux of dust species 1 out of layer [kg]
    real(r8) :: qin_dst2(lbc:ubc)                 ! flux of dust species 2 into layer [kg]
    real(r8) :: qout_dst2(lbc:ubc)                ! flux of dust species 2 out of layer [kg]
    real(r8) :: qin_dst3(lbc:ubc)                 ! flux of dust species 3 into layer [kg]
    real(r8) :: qout_dst3(lbc:ubc)                ! flux of dust species 3 out of layer [kg]
    real(r8) :: qin_dst4(lbc:ubc)                 ! flux of dust species 4 into layer [kg]
    real(r8) :: qout_dst4(lbc:ubc)                ! flux of dust species 4 out of layer [kg]
    real(r8) :: mss_liqice                        ! mass of liquid+ice in a layer
 
!-----------------------------------------------------------------------

    ! Assign local pointers to derived subtype components (column-level)

    snl              => clm3%g%l%c%cps%snl
    do_capsnow       => clm3%g%l%c%cps%do_capsnow
    qflx_snomelt     => clm3%g%l%c%cwf%qflx_snomelt
    qflx_rain_grnd   => clm3%g%l%c%cwf%pwf_a%qflx_rain_grnd
    qflx_sub_snow    => clm3%g%l%c%cwf%pwf_a%qflx_sub_snow
    qflx_evap_grnd   => clm3%g%l%c%cwf%pwf_a%qflx_evap_grnd
    qflx_dew_snow    => clm3%g%l%c%cwf%pwf_a%qflx_dew_snow
    qflx_dew_grnd    => clm3%g%l%c%cwf%pwf_a%qflx_dew_grnd
    qflx_top_soil    => clm3%g%l%c%cwf%qflx_top_soil
    dz               => clm3%g%l%c%cps%dz
    h2osoi_ice       => clm3%g%l%c%cws%h2osoi_ice
    h2osoi_liq       => clm3%g%l%c%cws%h2osoi_liq
    cgridcell        => clm3%g%l%c%gridcell
    mss_bcphi        => clm3%g%l%c%cps%mss_bcphi
    mss_bcpho        => clm3%g%l%c%cps%mss_bcpho
    mss_ocphi        => clm3%g%l%c%cps%mss_ocphi
    mss_ocpho        => clm3%g%l%c%cps%mss_ocpho
    mss_dst1         => clm3%g%l%c%cps%mss_dst1
    mss_dst2         => clm3%g%l%c%cps%mss_dst2
    mss_dst3         => clm3%g%l%c%cps%mss_dst3
    mss_dst4         => clm3%g%l%c%cps%mss_dst4
    flx_bc_dep       => clm3%g%l%c%cwf%flx_bc_dep
    flx_bc_dep_wet   => clm3%g%l%c%cwf%flx_bc_dep_wet
    flx_bc_dep_dry   => clm3%g%l%c%cwf%flx_bc_dep_dry
    flx_bc_dep_phi   => clm3%g%l%c%cwf%flx_bc_dep_phi
    flx_bc_dep_pho   => clm3%g%l%c%cwf%flx_bc_dep_pho
    flx_oc_dep       => clm3%g%l%c%cwf%flx_oc_dep
    flx_oc_dep_wet   => clm3%g%l%c%cwf%flx_oc_dep_wet
    flx_oc_dep_dry   => clm3%g%l%c%cwf%flx_oc_dep_dry
    flx_oc_dep_phi   => clm3%g%l%c%cwf%flx_oc_dep_phi
    flx_oc_dep_pho   => clm3%g%l%c%cwf%flx_oc_dep_pho
    flx_dst_dep      => clm3%g%l%c%cwf%flx_dst_dep
    flx_dst_dep_wet1 => clm3%g%l%c%cwf%flx_dst_dep_wet1
    flx_dst_dep_dry1 => clm3%g%l%c%cwf%flx_dst_dep_dry1
    flx_dst_dep_wet2 => clm3%g%l%c%cwf%flx_dst_dep_wet2
    flx_dst_dep_dry2 => clm3%g%l%c%cwf%flx_dst_dep_dry2
    flx_dst_dep_wet3 => clm3%g%l%c%cwf%flx_dst_dep_wet3
    flx_dst_dep_dry3 => clm3%g%l%c%cwf%flx_dst_dep_dry3
    flx_dst_dep_wet4 => clm3%g%l%c%cwf%flx_dst_dep_wet4
    flx_dst_dep_dry4 => clm3%g%l%c%cwf%flx_dst_dep_dry4
    forc_aer         => clm_a2l%forc_aer

    ! Determine model time step

    dtime = get_step_size()

    ! Renew the mass of ice lens (h2osoi_ice) and liquid (h2osoi_liq) in the
    ! surface snow layer resulting from sublimation (frost) / evaporation (condense)

    do fc = 1,num_snowc
       c = filter_snowc(fc)
       if (do_capsnow(c)) then
          wgdif = h2osoi_ice(c,snl(c)+1) - qflx_sub_snow(c)*dtime
          h2osoi_ice(c,snl(c)+1) = wgdif
          if (wgdif < 0._r8) then
             h2osoi_ice(c,snl(c)+1) = 0._r8
             h2osoi_liq(c,snl(c)+1) = h2osoi_liq(c,snl(c)+1) + wgdif
          end if
          h2osoi_liq(c,snl(c)+1) = h2osoi_liq(c,snl(c)+1) - qflx_evap_grnd(c) * dtime
       else
          wgdif = h2osoi_ice(c,snl(c)+1) + (qflx_dew_snow(c) - qflx_sub_snow(c)) * dtime
          h2osoi_ice(c,snl(c)+1) = wgdif
          if (wgdif < 0._r8) then
             h2osoi_ice(c,snl(c)+1) = 0._r8
             h2osoi_liq(c,snl(c)+1) = h2osoi_liq(c,snl(c)+1) + wgdif
          end if
          h2osoi_liq(c,snl(c)+1) = h2osoi_liq(c,snl(c)+1) +  &
               (qflx_rain_grnd(c) + qflx_dew_grnd(c) - qflx_evap_grnd(c)) * dtime
       end if
       h2osoi_liq(c,snl(c)+1) = max(0._r8, h2osoi_liq(c,snl(c)+1))
    end do

    ! Porosity and partial volume

    do j = -nlevsno+1, 0
       do fc = 1, num_snowc
          c = filter_snowc(fc)
          if (j >= snl(c)+1) then
             vol_ice(c,j) = min(1._r8, h2osoi_ice(c,j)/(dz(c,j)*denice))
             eff_porosity(c,j) = 1._r8 - vol_ice(c,j)
             vol_liq(c,j) = min(eff_porosity(c,j),h2osoi_liq(c,j)/(dz(c,j)*denh2o))
          end if
       end do
    end do

    ! Capillary forces within snow are usually two or more orders of magnitude
    ! less than those of gravity. Only gravity terms are considered.
    ! the genernal expression for water flow is "K * ss**3", however,
    ! no effective parameterization for "K".  Thus, a very simple consideration
    ! (not physically based) is introduced:
    ! when the liquid water of layer exceeds the layer's holding
    ! capacity, the excess meltwater adds to the underlying neighbor layer.

    ! Also compute aerosol fluxes through snowpack in this loop:
    ! 1) compute aerosol mass in each layer
    ! 2) add aerosol mass flux from above layer to mass of this layer
    ! 3) qout_xxx is mass flux of aerosol species xxx out bottom of
    !    layer in water flow, proportional to (current) concentration
    !    of aerosol in layer multiplied by a scavenging ratio.
    ! 4) update mass of aerosol in top layer, accordingly
    ! 5) update mass concentration of aerosol accordingly

    qin(:) = 0._r8
    qin_bc_phi(:) = 0._r8
    qin_bc_pho(:) = 0._r8
    qin_oc_phi(:) = 0._r8
    qin_oc_pho(:) = 0._r8
    qin_dst1(:)   = 0._r8
    qin_dst2(:)   = 0._r8
    qin_dst3(:)   = 0._r8
    qin_dst4(:)   = 0._r8

    do j = -nlevsno+1, 0
       do fc = 1, num_snowc
          c = filter_snowc(fc)
          if (j >= snl(c)+1) then
             h2osoi_liq(c,j) = h2osoi_liq(c,j) + qin(c)
             
             mss_bcphi(c,j) = mss_bcphi(c,j) + qin_bc_phi(c)
             mss_bcpho(c,j) = mss_bcpho(c,j) + qin_bc_pho(c)
             mss_ocphi(c,j) = mss_ocphi(c,j) + qin_oc_phi(c)
             mss_ocpho(c,j) = mss_ocpho(c,j) + qin_oc_pho(c)
             mss_dst1(c,j)  = mss_dst1(c,j) + qin_dst1(c)
             mss_dst2(c,j)  = mss_dst2(c,j) + qin_dst2(c)
             mss_dst3(c,j)  = mss_dst3(c,j) + qin_dst3(c)
             mss_dst4(c,j)  = mss_dst4(c,j) + qin_dst4(c)

             if (j <= -1) then
                ! No runoff over snow surface, just ponding on surface
                if (eff_porosity(c,j) < wimp .OR. eff_porosity(c,j+1) < wimp) then
                   qout(c) = 0._r8
                else
                   qout(c) = max(0._r8,(vol_liq(c,j)-ssi*eff_porosity(c,j))*dz(c,j))
                   qout(c) = min(qout(c),(1._r8-vol_ice(c,j+1)-vol_liq(c,j+1))*dz(c,j+1))
                end if
             else
                qout(c) = max(0._r8,(vol_liq(c,j) - ssi*eff_porosity(c,j))*dz(c,j))
             end if
             qout(c) = qout(c)*1000._r8
             h2osoi_liq(c,j) = h2osoi_liq(c,j) - qout(c)
             qin(c) = qout(c)

             ! mass of ice+water: in extremely rare circumstances, this can
             ! be zero, even though there is a snow layer defined. In
             ! this case, set the mass to a very small value to
             ! prevent division by zero.
             mss_liqice = h2osoi_liq(c,j)+h2osoi_ice(c,j)
             if (mss_liqice < 1E-30_r8) then
                mss_liqice = 1E-30_r8
             endif
         
             ! BCPHI:
             ! 1. flux with meltwater:
             qout_bc_phi(c) = qout(c)*scvng_fct_mlt_bcphi*(mss_bcphi(c,j)/mss_liqice)
             if (qout_bc_phi(c) > mss_bcphi(c,j)) then
                qout_bc_phi(c) = mss_bcphi(c,j)
             endif
             mss_bcphi(c,j) = mss_bcphi(c,j) - qout_bc_phi(c)
             qin_bc_phi(c) = qout_bc_phi(c)

             ! BCPHO:
             ! 1. flux with meltwater:
             qout_bc_pho(c) = qout(c)*scvng_fct_mlt_bcpho*(mss_bcpho(c,j)/mss_liqice)
             if (qout_bc_pho(c) > mss_bcpho(c,j)) then
                qout_bc_pho(c) = mss_bcpho(c,j)
             endif
             mss_bcpho(c,j) = mss_bcpho(c,j) - qout_bc_pho(c)
             qin_bc_pho(c) = qout_bc_pho(c)

             ! OCPHI:
             ! 1. flux with meltwater:
             qout_oc_phi(c) = qout(c)*scvng_fct_mlt_ocphi*(mss_ocphi(c,j)/mss_liqice)
             if (qout_oc_phi(c) > mss_ocphi(c,j)) then
                qout_oc_phi(c) = mss_ocphi(c,j)
             endif
             mss_ocphi(c,j) = mss_ocphi(c,j) - qout_oc_phi(c)
             qin_oc_phi(c) = qout_oc_phi(c)

             ! OCPHO:
             ! 1. flux with meltwater:
             qout_oc_pho(c) = qout(c)*scvng_fct_mlt_ocpho*(mss_ocpho(c,j)/mss_liqice)
             if (qout_oc_pho(c) > mss_ocpho(c,j)) then
                qout_oc_pho(c) = mss_ocpho(c,j)
             endif
             mss_ocpho(c,j) = mss_ocpho(c,j) - qout_oc_pho(c)
             qin_oc_pho(c) = qout_oc_pho(c)

             ! DUST 1:
             ! 1. flux with meltwater:
             qout_dst1(c) = qout(c)*scvng_fct_mlt_dst1*(mss_dst1(c,j)/mss_liqice)
             if (qout_dst1(c) > mss_dst1(c,j)) then
                qout_dst1(c) = mss_dst1(c,j)
             endif
             mss_dst1(c,j) = mss_dst1(c,j) - qout_dst1(c)
             qin_dst1(c) = qout_dst1(c)

             ! DUST 2:
             ! 1. flux with meltwater:
             qout_dst2(c) = qout(c)*scvng_fct_mlt_dst2*(mss_dst2(c,j)/mss_liqice)
             if (qout_dst2(c) > mss_dst2(c,j)) then
                qout_dst2(c) = mss_dst2(c,j)
             endif
             mss_dst2(c,j) = mss_dst2(c,j) - qout_dst2(c)
             qin_dst2(c) = qout_dst2(c)

             ! DUST 3:
             ! 1. flux with meltwater:
             qout_dst3(c) = qout(c)*scvng_fct_mlt_dst3*(mss_dst3(c,j)/mss_liqice)
             if (qout_dst3(c) > mss_dst3(c,j)) then
                qout_dst3(c) = mss_dst3(c,j)
             endif
             mss_dst3(c,j) = mss_dst3(c,j) - qout_dst3(c)
             qin_dst3(c) = qout_dst3(c)

             ! DUST 4:
             ! 1. flux with meltwater:
             qout_dst4(c) = qout(c)*scvng_fct_mlt_dst4*(mss_dst4(c,j)/mss_liqice)
             if (qout_dst4(c) > mss_dst4(c,j)) then
                qout_dst4(c) = mss_dst4(c,j)
             endif
             mss_dst4(c,j) = mss_dst4(c,j) - qout_dst4(c)
             qin_dst4(c) = qout_dst4(c)
             
          end if
       end do
    end do

    ! Adjust layer thickness for any water+ice content changes in excess of previous 
    ! layer thickness. Strictly speaking, only necessary for top snow layer, but doing
    ! it for all snow layers will catch problems with older initial files.
    ! Layer interfaces (zi) and node depths (z) do not need adjustment here because they
    ! are adjusted in CombineSnowLayers and are not used up to that point.

    do j = -nlevsno+1, 0
       do fc = 1, num_snowc
          c = filter_snowc(fc)
          if (j >= snl(c)+1) then
              dz(c,j) = max(dz(c,j),h2osoi_liq(c,j)/denh2o + h2osoi_ice(c,j)/denice)
          end if
       end do
    end do
 
    do fc = 1, num_snowc
       c = filter_snowc(fc)
       ! Qout from snow bottom
       qflx_top_soil(c) = qout(c) / dtime
    end do

    do fc = 1, num_nosnowc
       c = filter_nosnowc(fc)
       qflx_top_soil(c) = qflx_rain_grnd(c) + qflx_snomelt(c)
    end do

    
    !  set aerosol deposition fluxes from forcing array
    !  The forcing array is either set from an external file 
    !  or from fluxes received from the atmosphere model
    do c = lbc,ubc
       g = cgridcell(c)
       
       flx_bc_dep_dry(c)   = forc_aer(g,1) + forc_aer(g,2)
       flx_bc_dep_wet(c)   = forc_aer(g,3)
       flx_bc_dep_phi(c)   = forc_aer(g,1) + forc_aer(g,3)
       flx_bc_dep_pho(c)   = forc_aer(g,2)
       flx_bc_dep(c)       = forc_aer(g,1) + forc_aer(g,2) + forc_aer(g,3)
       
       flx_oc_dep_dry(c)   = forc_aer(g,4) + forc_aer(g,5)
       flx_oc_dep_wet(c)   = forc_aer(g,6)
       flx_oc_dep_phi(c)   = forc_aer(g,4) + forc_aer(g,6)
       flx_oc_dep_pho(c)   = forc_aer(g,5)
       flx_oc_dep(c)       = forc_aer(g,4) + forc_aer(g,5) + forc_aer(g,6)
       
       flx_dst_dep_wet1(c) = forc_aer(g,7)
       flx_dst_dep_dry1(c) = forc_aer(g,8)
       flx_dst_dep_wet2(c) = forc_aer(g,9)
       flx_dst_dep_dry2(c) = forc_aer(g,10)
       flx_dst_dep_wet3(c) = forc_aer(g,11)
       flx_dst_dep_dry3(c) = forc_aer(g,12)
       flx_dst_dep_wet4(c) = forc_aer(g,13)
       flx_dst_dep_dry4(c) = forc_aer(g,14)
       flx_dst_dep(c)      = forc_aer(g,7) + forc_aer(g,8) + forc_aer(g,9) + &
                             forc_aer(g,10) + forc_aer(g,11) + forc_aer(g,12) + &
                             forc_aer(g,13) + forc_aer(g,14)
    
    end do

    ! aerosol deposition fluxes into top layer
    ! This is done after the inter-layer fluxes so that some aerosol
    ! is in the top layer after deposition, and is not immediately
    ! washed out before radiative calculations are done
    do fc = 1, num_snowc
       c = filter_snowc(fc)
       mss_bcphi(c,snl(c)+1) = mss_bcphi(c,snl(c)+1) + (flx_bc_dep_phi(c)*dtime)
       mss_bcpho(c,snl(c)+1) = mss_bcpho(c,snl(c)+1) + (flx_bc_dep_pho(c)*dtime)
       mss_ocphi(c,snl(c)+1) = mss_ocphi(c,snl(c)+1) + (flx_oc_dep_phi(c)*dtime)
       mss_ocpho(c,snl(c)+1) = mss_ocpho(c,snl(c)+1) + (flx_oc_dep_pho(c)*dtime)
       
       mss_dst1(c,snl(c)+1) = mss_dst1(c,snl(c)+1) + (flx_dst_dep_dry1(c) + flx_dst_dep_wet1(c))*dtime
       mss_dst2(c,snl(c)+1) = mss_dst2(c,snl(c)+1) + (flx_dst_dep_dry2(c) + flx_dst_dep_wet2(c))*dtime
       mss_dst3(c,snl(c)+1) = mss_dst3(c,snl(c)+1) + (flx_dst_dep_dry3(c) + flx_dst_dep_wet3(c))*dtime
       mss_dst4(c,snl(c)+1) = mss_dst4(c,snl(c)+1) + (flx_dst_dep_dry4(c) + flx_dst_dep_wet4(c))*dtime
    end do

  end subroutine SnowWater

!-----------------------------------------------------------------------
!BOP
!
! !IROUTINE: SnowCompaction
!
! !INTERFACE:
  subroutine SnowCompaction(lbc, ubc, num_snowc, filter_snowc)
!
! !DESCRIPTION:
! Determine the change in snow layer thickness due to compaction and
! settling.
! Three metamorphisms of changing snow characteristics are implemented,
! i.e., destructive, overburden, and melt. The treatments of the former
! two are from SNTHERM.89 and SNTHERM.99 (1991, 1999). The contribution
! due to melt metamorphism is simply taken as a ratio of snow ice
! fraction after the melting versus before the melting.
!
! !USES:
    use clmtype
    use clm_time_manager, only : get_step_size
    use clm_varcon  , only : denice, denh2o, tfrz, istice_mec
!
! !ARGUMENTS:
    implicit none
    integer, intent(in) :: lbc, ubc                ! column bounds
    integer, intent(in) :: num_snowc               ! number of column snow points in column filter
    integer, intent(in) :: filter_snowc(ubc-lbc+1) ! column filter for snow points
!
! !CALLED FROM:
! subroutine Hydrology2 in module Hydrology2Mod
!
! !REVISION HISTORY:
! 15 September 1999: Yongjiu Dai; Initial code
! 15 December 1999:  Paul Houser and Jon Radakovich; F90 Revision
! 2/28/02, Peter Thornton: Migrated to new data structures
! 2/29/08, David Lawrence: Revised snow overburden to be include 0.5 weight of current layer
!
! !LOCAL VARIABLES:
!
! local pointers to implicit in scalars
!
    integer,  pointer :: snl(:)             !number of snow layers
!
! local pointers to implicit in arguments
!
    integer,  pointer :: imelt(:,:)        !flag for melting (=1), freezing (=2), Not=0
    real(r8), pointer :: frac_iceold(:,:)  !fraction of ice relative to the tot water
    real(r8), pointer :: t_soisno(:,:)     !soil temperature (Kelvin)
    real(r8), pointer :: h2osoi_ice(:,:)   !ice lens (kg/m2)
    real(r8), pointer :: h2osoi_liq(:,:)   !liquid water (kg/m2)
!
! local pointers to implicit inout arguments
!
    real(r8), pointer :: dz(:,:)           !layer depth (m)
!
!
! !OTHER LOCAL VARIABLES:
!EOP
!
    integer :: j, l, c, fc                   ! indices
    real(r8):: dtime                         ! land model time step (sec)
    real(r8), parameter :: c2 = 23.e-3_r8    ! [m3/kg]
    real(r8), parameter :: c3 = 2.777e-6_r8  ! [1/s]
    real(r8), parameter :: c4 = 0.04_r8      ! [1/K]
    real(r8), parameter :: c5 = 2.0_r8       !
    real(r8), parameter :: dm = 100.0_r8     ! Upper Limit on Destructive Metamorphism Compaction [kg/m3]
    real(r8), parameter :: eta0 = 9.e+5_r8   ! The Viscosity Coefficient Eta0 [kg-s/m2]
    real(r8) :: burden(lbc:ubc) ! pressure of overlying snow [kg/m2]
    real(r8) :: ddz1   ! Rate of settling of snowpack due to destructive metamorphism.
    real(r8) :: ddz2   ! Rate of compaction of snowpack due to overburden.
    real(r8) :: ddz3   ! Rate of compaction of snowpack due to melt [1/s]
    real(r8) :: dexpf  ! expf=exp(-c4*(273.15-t_soisno)).
    real(r8) :: fi     ! Fraction of ice relative to the total water content at current time step
    real(r8) :: td     ! t_soisno - tfrz [K]
    real(r8) :: pdzdtc ! Nodal rate of change in fractional-thickness due to compaction [fraction/s]
    real(r8) :: void   ! void (1 - vol_ice - vol_liq)
    real(r8) :: wx     ! water mass (ice+liquid) [kg/m2]
    real(r8) :: bi     ! partial density of ice [kg/m3]

    integer, pointer :: clandunit(:)       !landunit index for each column
    integer, pointer :: ltype(:)           !landunit type

!-----------------------------------------------------------------------

    ! Assign local pointers to derived subtypes (column-level)

    snl         => clm3%g%l%c%cps%snl
    dz          => clm3%g%l%c%cps%dz
    imelt       => clm3%g%l%c%cps%imelt
    frac_iceold => clm3%g%l%c%cps%frac_iceold
    t_soisno    => clm3%g%l%c%ces%t_soisno
    h2osoi_ice  => clm3%g%l%c%cws%h2osoi_ice
    h2osoi_liq  => clm3%g%l%c%cws%h2osoi_liq
    clandunit   => clm3%g%l%c%landunit
    ltype       => clm3%g%l%itype

    ! Get time step

    dtime = get_step_size()

    ! Begin calculation - note that the following column loops are only invoked if snl(c) < 0

    burden(:) = 0._r8

    do j = -nlevsno+1, 0
       do fc = 1, num_snowc
          c = filter_snowc(fc)
          if (j >= snl(c)+1) then

             wx = h2osoi_ice(c,j) + h2osoi_liq(c,j)
             void = 1._r8 - (h2osoi_ice(c,j)/denice + h2osoi_liq(c,j)/denh2o) / dz(c,j)

             ! If void is negative, then increase dz such that void = 0.
             ! This should be done for any landunit, but for now is done only for glacier_mec 1andunits.
             l = clandunit(c)
             if (ltype(l)==istice_mec .and. void < 0._r8) then
                dz(c,j) = h2osoi_ice(c,j)/denice + h2osoi_liq(c,j)/denh2o
                void = 0._r8
             endif

             ! Allow compaction only for non-saturated node and higher ice lens node.
             if (void > 0.001_r8 .and. h2osoi_ice(c,j) > .1_r8) then
                bi = h2osoi_ice(c,j) / dz(c,j)
                fi = h2osoi_ice(c,j) / wx
                td = tfrz-t_soisno(c,j)
                dexpf = exp(-c4*td)

                ! Settling as a result of destructive metamorphism

                ddz1 = -c3*dexpf
                if (bi > dm) ddz1 = ddz1*exp(-46.0e-3_r8*(bi-dm))

                ! Liquid water term

                if (h2osoi_liq(c,j) > 0.01_r8*dz(c,j)) ddz1=ddz1*c5

                ! Compaction due to overburden

                ddz2 = -(burden(c)+wx/2._r8)*exp(-0.08_r8*td - c2*bi)/eta0 

                ! Compaction occurring during melt

                if (imelt(c,j) == 1) then
                   ddz3 = - 1._r8/dtime * max(0._r8,(frac_iceold(c,j) - fi)/frac_iceold(c,j))
                else
                   ddz3 = 0._r8
                end if

                ! Time rate of fractional change in dz (units of s-1)

                pdzdtc = ddz1 + ddz2 + ddz3

                ! The change in dz due to compaction
                ! Limit compaction to no less than fully saturated layer thickness

                dz(c,j) = max(dz(c,j) * (1._r8+pdzdtc*dtime),h2osoi_ice(c,j)/denice &
                                         + h2osoi_liq(c,j)/denh2o)

             end if

             ! Pressure of overlying snow

             burden(c) = burden(c) + wx

          end if
       end do
    end do

  end subroutine SnowCompaction

!-----------------------------------------------------------------------
!BOP
!
! !IROUTINE: CombineSnowLayers
!
! !INTERFACE:
  subroutine CombineSnowLayers(lbc, ubc, num_snowc, filter_snowc)
!
! !DESCRIPTION:
! Combine snow layers that are less than a minimum thickness or mass
! If the snow element thickness or mass is less than a prescribed minimum,
! then it is combined with a neighboring element.  The subroutine
! clm\_combo.f90 then executes the combination of mass and energy.
!
! !USES:
    use clmtype
    use clm_varcon, only : istsoil, isturb
    use clm_varcon, only : istcrop
!
! !ARGUMENTS:
    implicit none
    integer, intent(in)    :: lbc, ubc                    ! column bounds
    integer, intent(inout) :: num_snowc                   ! number of column snow points in column filter
    integer, intent(inout) :: filter_snowc(ubc-lbc+1)     ! column filter for snow points
!
! !CALLED FROM:
! subroutine Hydrology2 in module Hydrology2Mod
!
! !REVISION HISTORY:
! 15 September 1999: Yongjiu Dai; Initial code
! 15 December 1999:  Paul Houser and Jon Radakovich; F90 Revision
! 2/28/02, Peter Thornton: Migrated to new data structures.
! 03/28/08, Mark Flanner: Added aerosol masses and snow grain radius
!
! !LOCAL VARIABLES:
!
! local pointers to implicit in arguments
!
    integer, pointer :: clandunit(:)       !landunit index for each column
    integer, pointer :: ltype(:)           !landunit type
!
! local pointers to implicit inout arguments
!
    integer , pointer :: snl(:)            !number of snow layers
    real(r8), pointer :: h2osno(:)         !snow water (mm H2O)
    real(r8), pointer :: snowdp(:)         !snow height (m)
    real(r8), pointer :: dz(:,:)           !layer depth (m)
    real(r8), pointer :: zi(:,:)           !interface level below a "z" level (m)
    real(r8), pointer :: t_soisno(:,:)     !soil temperature (Kelvin)
    real(r8), pointer :: h2osoi_ice(:,:)   !ice lens (kg/m2)
    real(r8), pointer :: h2osoi_liq(:,:)   !liquid water (kg/m2)
!
! local pointers to implicit out arguments
!
    real(r8), pointer :: z(:,:)          ! layer thickness (m)
    real(r8), pointer :: mss_bcphi(:,:)  ! hydrophilic BC mass in snow (col,lyr) [kg]
    real(r8), pointer :: mss_bcpho(:,:)  ! hydrophobic BC mass in snow (col,lyr) [kg]
    real(r8), pointer :: mss_ocphi(:,:)  ! hydrophilic OC mass in snow (col,lyr) [kg]
    real(r8), pointer :: mss_ocpho(:,:)  ! hydrophobic OC mass in snow (col,lyr) [kg]
    real(r8), pointer :: mss_dst1(:,:)   ! dust species 1 mass in snow (col,lyr) [kg]
    real(r8), pointer :: mss_dst2(:,:)   ! dust species 2 mass in snow (col,lyr) [kg]
    real(r8), pointer :: mss_dst3(:,:)   ! dust species 3 mass in snow (col,lyr) [kg]
    real(r8), pointer :: mss_dst4(:,:)   ! dust species 4 mass in snow (col,lyr) [kg]
    real(r8), pointer :: snw_rds(:,:)    ! effective snow grain radius (col,lyr) [microns, m^-6]
!
!
! !OTHER LOCAL VARIABLES:
!EOP
!
    integer :: c, fc                 ! column indices
    integer :: i,k                   ! loop indices
    integer :: j,l                   ! node indices
    integer :: msn_old(lbc:ubc)      ! number of top snow layer
    integer :: mssi(lbc:ubc)         ! node index
    integer :: neibor                ! adjacent node selected for combination
    real(r8):: zwice(lbc:ubc)        ! total ice mass in snow
    real(r8):: zwliq (lbc:ubc)       ! total liquid water in snow
    real(r8):: dzmin(5)              ! minimum of top snow layer

    data dzmin /0.010_r8, 0.015_r8, 0.025_r8, 0.055_r8, 0.115_r8/
!-----------------------------------------------------------------------

    ! Assign local pointers to derived subtypes (landunit-level)

    ltype      => clm3%g%l%itype

    ! Assign local pointers to derived subtypes (column-level)

    clandunit  => clm3%g%l%c%landunit
    snl        => clm3%g%l%c%cps%snl
    snowdp     => clm3%g%l%c%cps%snowdp
    h2osno     => clm3%g%l%c%cws%h2osno
    dz         => clm3%g%l%c%cps%dz
    zi         => clm3%g%l%c%cps%zi
    z          => clm3%g%l%c%cps%z
    t_soisno   => clm3%g%l%c%ces%t_soisno
    h2osoi_ice => clm3%g%l%c%cws%h2osoi_ice
    h2osoi_liq => clm3%g%l%c%cws%h2osoi_liq
    mss_bcphi  => clm3%g%l%c%cps%mss_bcphi
    mss_bcpho  => clm3%g%l%c%cps%mss_bcpho
    mss_ocphi  => clm3%g%l%c%cps%mss_ocphi
    mss_ocpho  => clm3%g%l%c%cps%mss_ocpho
    mss_dst1   => clm3%g%l%c%cps%mss_dst1
    mss_dst2   => clm3%g%l%c%cps%mss_dst2
    mss_dst3   => clm3%g%l%c%cps%mss_dst3
    mss_dst4   => clm3%g%l%c%cps%mss_dst4
    snw_rds    => clm3%g%l%c%cps%snw_rds


    ! Check the mass of ice lens of snow, when the total is less than a small value,
    ! combine it with the underlying neighbor.

    do fc = 1, num_snowc
       c = filter_snowc(fc)
       msn_old(c) = snl(c)
    end do

    ! The following loop is NOT VECTORIZED

    do fc = 1, num_snowc
       c = filter_snowc(fc)
       l = clandunit(c)
       do j = msn_old(c)+1,0
          if (h2osoi_ice(c,j) <= .1_r8) then
             if (ltype(l) == istsoil .or. ltype(l)==isturb .or. ltype(l) == istcrop) then
                h2osoi_liq(c,j+1) = h2osoi_liq(c,j+1) + h2osoi_liq(c,j)
                h2osoi_ice(c,j+1) = h2osoi_ice(c,j+1) + h2osoi_ice(c,j)

                if (j /= 0) dz(c,j+1) = dz(c,j+1) + dz(c,j)
                
                ! NOTE: Temperature, and similarly snw_rds, of the
                ! underlying snow layer are NOT adjusted in this case. 
                ! Because the layer being eliminated has a small mass, 
                ! this should not make a large difference, but it 
                ! would be more thorough to do so.
                if (j /= 0) then
                   mss_bcphi(c,j+1) = mss_bcphi(c,j+1) + mss_bcphi(c,j)
                   mss_bcpho(c,j+1) = mss_bcpho(c,j+1) + mss_bcpho(c,j)
                   mss_ocphi(c,j+1) = mss_ocphi(c,j+1) + mss_ocphi(c,j)
                   mss_ocpho(c,j+1) = mss_ocpho(c,j+1) + mss_ocpho(c,j)
                   mss_dst1(c,j+1)  = mss_dst1(c,j+1) + mss_dst1(c,j)
                   mss_dst2(c,j+1)  = mss_dst2(c,j+1) + mss_dst2(c,j)
                   mss_dst3(c,j+1)  = mss_dst3(c,j+1) + mss_dst3(c,j)
                   mss_dst4(c,j+1)  = mss_dst4(c,j+1) + mss_dst4(c,j)
                end if

             else if (ltype(l) /= istsoil .and. ltype(l) /= isturb .and. ltype(l) /= istcrop .and. j /= 0) then
                h2osoi_liq(c,j+1) = h2osoi_liq(c,j+1) + h2osoi_liq(c,j)
                h2osoi_ice(c,j+1) = h2osoi_ice(c,j+1) + h2osoi_ice(c,j)
                dz(c,j+1) = dz(c,j+1) + dz(c,j)
                
                mss_bcphi(c,j+1) = mss_bcphi(c,j+1) + mss_bcphi(c,j)
                mss_bcpho(c,j+1) = mss_bcpho(c,j+1) + mss_bcpho(c,j)
                mss_ocphi(c,j+1) = mss_ocphi(c,j+1) + mss_ocphi(c,j)
                mss_ocpho(c,j+1) = mss_ocpho(c,j+1) + mss_ocpho(c,j)
                mss_dst1(c,j+1)  = mss_dst1(c,j+1) + mss_dst1(c,j)
                mss_dst2(c,j+1)  = mss_dst2(c,j+1) + mss_dst2(c,j)
                mss_dst3(c,j+1)  = mss_dst3(c,j+1) + mss_dst3(c,j)
                mss_dst4(c,j+1)  = mss_dst4(c,j+1) + mss_dst4(c,j)

             end if

             ! shift all elements above this down one.
             if (j > snl(c)+1 .and. snl(c) < -1) then
                do i = j, snl(c)+2, -1
                   t_soisno(c,i)   = t_soisno(c,i-1)
                   h2osoi_liq(c,i) = h2osoi_liq(c,i-1)
                   h2osoi_ice(c,i) = h2osoi_ice(c,i-1)
                   
                   mss_bcphi(c,i)   = mss_bcphi(c,i-1)
                   mss_bcpho(c,i)   = mss_bcpho(c,i-1)
                   mss_ocphi(c,i)   = mss_ocphi(c,i-1)
                   mss_ocpho(c,i)   = mss_ocpho(c,i-1)
                   mss_dst1(c,i)    = mss_dst1(c,i-1)
                   mss_dst2(c,i)    = mss_dst2(c,i-1)
                   mss_dst3(c,i)    = mss_dst3(c,i-1)
                   mss_dst4(c,i)    = mss_dst4(c,i-1)
                   snw_rds(c,i)     = snw_rds(c,i-1)

                   dz(c,i)         = dz(c,i-1)
                end do
             end if
             snl(c) = snl(c) + 1
          end if
       end do
    end do

    do fc = 1, num_snowc
       c = filter_snowc(fc)
       h2osno(c) = 0._r8
       snowdp(c) = 0._r8
       zwice(c)  = 0._r8
       zwliq(c)  = 0._r8
    end do

    do j = -nlevsno+1,0
       do fc = 1, num_snowc
          c = filter_snowc(fc)
          if (j >= snl(c)+1) then
             h2osno(c) = h2osno(c) + h2osoi_ice(c,j) + h2osoi_liq(c,j)
             snowdp(c) = snowdp(c) + dz(c,j)
             zwice(c)  = zwice(c) + h2osoi_ice(c,j)
             zwliq(c)  = zwliq(c) + h2osoi_liq(c,j)
          end if
       end do
    end do

    ! Check the snow depth - all snow gone
    ! The liquid water assumes ponding on soil surface.

    do fc = 1, num_snowc
       c = filter_snowc(fc)
       l = clandunit(c)
       if (snowdp(c) < 0.01_r8 .and. snowdp(c) > 0._r8) then
          snl(c) = 0
          h2osno(c) = zwice(c)

          mss_bcphi(c,:) = 0._r8
          mss_bcpho(c,:) = 0._r8
          mss_ocphi(c,:) = 0._r8
          mss_ocpho(c,:) = 0._r8
          mss_dst1(c,:)  = 0._r8
          mss_dst2(c,:)  = 0._r8
          mss_dst3(c,:)  = 0._r8
          mss_dst4(c,:)  = 0._r8

          if (h2osno(c) <= 0._r8) snowdp(c) = 0._r8
          if (ltype(l) == istsoil .or. ltype(l) == isturb .or. ltype(l) == istcrop) then
             h2osoi_liq(c,1) = h2osoi_liq(c,1) + zwliq(c)
          end if
       end if
    end do

    ! Check the snow depth - snow layers combined
    ! The following loop IS NOT VECTORIZED

    do fc = 1, num_snowc
       c = filter_snowc(fc)

       ! Two or more layers

       if (snl(c) < -1) then

          msn_old(c) = snl(c)
          mssi(c) = 1

          do i = msn_old(c)+1,0
             if (dz(c,i) < dzmin(mssi(c))) then

                if (i == snl(c)+1) then
                   ! If top node is removed, combine with bottom neighbor.
                   neibor = i + 1
                else if (i == 0) then
                   ! If the bottom neighbor is not snow, combine with the top neighbor.
                   neibor = i - 1
                else
                   ! If none of the above special cases apply, combine with the thinnest neighbor
                   neibor = i + 1
                   if ((dz(c,i-1)+dz(c,i)) < (dz(c,i+1)+dz(c,i))) neibor = i-1
                end if

                ! Node l and j are combined and stored as node j.
                if (neibor > i) then
                   j = neibor
                   l = i
                else
                   j = i
                   l = neibor
                end if

                ! this should be included in 'Combo' for consistency,
                ! but functionally it is the same to do it here
                mss_bcphi(c,j)=mss_bcphi(c,j)+mss_bcphi(c,l)
                mss_bcpho(c,j)=mss_bcpho(c,j)+mss_bcpho(c,l)
                mss_ocphi(c,j)=mss_ocphi(c,j)+mss_ocphi(c,l)
                mss_ocpho(c,j)=mss_ocpho(c,j)+mss_ocpho(c,l)
                mss_dst1(c,j)=mss_dst1(c,j)+mss_dst1(c,l)
                mss_dst2(c,j)=mss_dst2(c,j)+mss_dst2(c,l)
                mss_dst3(c,j)=mss_dst3(c,j)+mss_dst3(c,l)
                mss_dst4(c,j)=mss_dst4(c,j)+mss_dst4(c,l)
                ! mass-weighted combination of effective grain size:
                snw_rds(c,j) = (snw_rds(c,j)*(h2osoi_liq(c,j)+h2osoi_ice(c,j)) + &
                               snw_rds(c,l)*(h2osoi_liq(c,l)+h2osoi_ice(c,l))) / &
                               (h2osoi_liq(c,j)+h2osoi_ice(c,j)+h2osoi_liq(c,l)+h2osoi_ice(c,l))

                call Combo (dz(c,j), h2osoi_liq(c,j), h2osoi_ice(c,j), &
                   t_soisno(c,j), dz(c,l), h2osoi_liq(c,l), h2osoi_ice(c,l), t_soisno(c,l) )

                ! Now shift all elements above this down one.
                if (j-1 > snl(c)+1) then
                   do k = j-1, snl(c)+2, -1
                      t_soisno(c,k) = t_soisno(c,k-1)
                      h2osoi_ice(c,k) = h2osoi_ice(c,k-1)
                      h2osoi_liq(c,k) = h2osoi_liq(c,k-1)

                      mss_bcphi(c,k) = mss_bcphi(c,k-1)
                      mss_bcpho(c,k) = mss_bcpho(c,k-1)
                      mss_ocphi(c,k) = mss_ocphi(c,k-1)
                      mss_ocpho(c,k) = mss_ocpho(c,k-1)
                      mss_dst1(c,k)  = mss_dst1(c,k-1)
                      mss_dst2(c,k)  = mss_dst2(c,k-1)
                      mss_dst3(c,k)  = mss_dst3(c,k-1)
                      mss_dst4(c,k)  = mss_dst4(c,k-1)
                      snw_rds(c,k)   = snw_rds(c,k-1)

                      dz(c,k) = dz(c,k-1)
                   end do
                end if

                ! Decrease the number of snow layers
                snl(c) = snl(c) + 1
                if (snl(c) >= -1) EXIT

             else

                ! The layer thickness is greater than the prescribed minimum value
                mssi(c) = mssi(c) + 1

             end if
          end do

       end if

    end do

    ! Reset the node depth and the depth of layer interface

    do j = 0, -nlevsno+1, -1
       do fc = 1, num_snowc
          c = filter_snowc(fc)
          if (j >= snl(c) + 1) then
             z(c,j) = zi(c,j) - 0.5_r8*dz(c,j)
             zi(c,j-1) = zi(c,j) - dz(c,j)
          end if
       end do
    end do

  end subroutine CombineSnowLayers

!-----------------------------------------------------------------------
!BOP
!
! !IROUTINE: DivideSnowLayers
!
! !INTERFACE:
  subroutine DivideSnowLayers(lbc, ubc, num_snowc, filter_snowc)
!
! !DESCRIPTION:
! Subdivides snow layers if they exceed their prescribed maximum thickness.
!
! !USES:
    use clmtype
    use clm_varcon,  only : tfrz 
!
! !ARGUMENTS:
    implicit none
    integer, intent(in)    :: lbc, ubc                    ! column bounds
    integer, intent(inout) :: num_snowc                   ! number of column snow points in column filter
    integer, intent(inout) :: filter_snowc(ubc-lbc+1)     ! column filter for snow points
!
! !CALLED FROM:
! subroutine Hydrology2 in module Hydrology2Mod
!
! !REVISION HISTORY:
! 15 September 1999: Yongjiu Dai; Initial code
! 15 December 1999:  Paul Houser and Jon Radakovich; F90 Revision
! 2/28/02, Peter Thornton: Migrated to new data structures.
! 2/29/08, David Lawrence: Snowpack T profile maintained during layer splitting
! 03/28/08, Mark Flanner: Added aerosol masses and snow grain radius
!
! !LOCAL VARIABLES:
!
! local pointers to implicit inout arguments
!
    integer , pointer :: snl(:)            !number of snow layers
    real(r8), pointer :: dz(:,:)           !layer depth (m)
    real(r8), pointer :: zi(:,:)           !interface level below a "z" level (m)
    real(r8), pointer :: t_soisno(:,:)     !soil temperature (Kelvin)
    real(r8), pointer :: h2osoi_ice(:,:)   !ice lens (kg/m2)
    real(r8), pointer :: h2osoi_liq(:,:)   !liquid water (kg/m2)
!
! local pointers to implicit out arguments
!
    real(r8), pointer :: z(:,:)          ! layer thickness (m)
    real(r8), pointer :: mss_bcphi(:,:)  ! hydrophilic BC mass in snow (col,lyr) [kg]
    real(r8), pointer :: mss_bcpho(:,:)  ! hydrophobic BC mass in snow (col,lyr) [kg]
    real(r8), pointer :: mss_ocphi(:,:)  ! hydrophilic OC mass in snow (col,lyr) [kg]
    real(r8), pointer :: mss_ocpho(:,:)  ! hydrophobic OC mass in snow (col,lyr) [kg]
    real(r8), pointer :: mss_dst1(:,:)   ! dust species 1 mass in snow (col,lyr) [kg]
    real(r8), pointer :: mss_dst2(:,:)   ! dust species 2 mass in snow (col,lyr) [kg]
    real(r8), pointer :: mss_dst3(:,:)   ! dust species 3 mass in snow (col,lyr) [kg]
    real(r8), pointer :: mss_dst4(:,:)   ! dust species 4 mass in snow (col,lyr) [kg]
    real(r8), pointer :: snw_rds(:,:)    ! effective snow grain radius (col,lyr) [microns, m^-6]
!
!
! !OTHER LOCAL VARIABLES:
!EOP
!
    integer  :: j, c, fc               ! indices
    real(r8) :: drr                    ! thickness of the combined [m]
    integer  :: msno                   ! number of snow layer 1 (top) to msno (bottom)
    real(r8) :: dzsno(lbc:ubc,nlevsno) ! Snow layer thickness [m]
    real(r8) :: swice(lbc:ubc,nlevsno) ! Partial volume of ice [m3/m3]
    real(r8) :: swliq(lbc:ubc,nlevsno) ! Partial volume of liquid water [m3/m3]
    real(r8) :: tsno(lbc:ubc ,nlevsno) ! Nodel temperature [K]
    real(r8) :: zwice                  ! temporary
    real(r8) :: zwliq                  ! temporary
    real(r8) :: propor                 ! temporary
    real(r8) :: dtdz                   ! temporary

    ! temporary variables mimicking the structure of other layer division variables
    real(r8) :: mbc_phi(lbc:ubc,nlevsno) ! mass of BC in each snow layer
    real(r8) :: zmbc_phi                 ! temporary
    real(r8) :: mbc_pho(lbc:ubc,nlevsno) ! mass of BC in each snow layer
    real(r8) :: zmbc_pho                 ! temporary
    real(r8) :: moc_phi(lbc:ubc,nlevsno) ! mass of OC in each snow layer
    real(r8) :: zmoc_phi                 ! temporary
    real(r8) :: moc_pho(lbc:ubc,nlevsno) ! mass of OC in each snow layer
    real(r8) :: zmoc_pho                 ! temporary
    real(r8) :: mdst1(lbc:ubc,nlevsno)   ! mass of dust 1 in each snow layer
    real(r8) :: zmdst1                   ! temporary
    real(r8) :: mdst2(lbc:ubc,nlevsno)   ! mass of dust 2 in each snow layer
    real(r8) :: zmdst2                   ! temporary
    real(r8) :: mdst3(lbc:ubc,nlevsno)   ! mass of dust 3 in each snow layer
    real(r8) :: zmdst3                   ! temporary
    real(r8) :: mdst4(lbc:ubc,nlevsno)   ! mass of dust 4 in each snow layer
    real(r8) :: zmdst4                   ! temporary
    real(r8) :: rds(lbc:ubc,nlevsno)

!-----------------------------------------------------------------------

    ! Assign local pointers to derived subtype components (column-level)

    snl        => clm3%g%l%c%cps%snl
    dz         => clm3%g%l%c%cps%dz
    zi         => clm3%g%l%c%cps%zi
    z          => clm3%g%l%c%cps%z
    t_soisno   => clm3%g%l%c%ces%t_soisno
    h2osoi_ice => clm3%g%l%c%cws%h2osoi_ice
    h2osoi_liq => clm3%g%l%c%cws%h2osoi_liq
    mss_bcphi  => clm3%g%l%c%cps%mss_bcphi
    mss_bcpho  => clm3%g%l%c%cps%mss_bcpho
    mss_ocphi  => clm3%g%l%c%cps%mss_ocphi
    mss_ocpho  => clm3%g%l%c%cps%mss_ocpho
    mss_dst1   => clm3%g%l%c%cps%mss_dst1
    mss_dst2   => clm3%g%l%c%cps%mss_dst2
    mss_dst3   => clm3%g%l%c%cps%mss_dst3
    mss_dst4   => clm3%g%l%c%cps%mss_dst4
    snw_rds    => clm3%g%l%c%cps%snw_rds
    

    ! Begin calculation - note that the following column loops are only invoked
    ! for snow-covered columns

    do j = 1,nlevsno
       do fc = 1, num_snowc
          c = filter_snowc(fc)
          if (j <= abs(snl(c))) then
             dzsno(c,j) = dz(c,j+snl(c))
             swice(c,j) = h2osoi_ice(c,j+snl(c))
             swliq(c,j) = h2osoi_liq(c,j+snl(c))
             tsno(c,j)  = t_soisno(c,j+snl(c))

             mbc_phi(c,j) = mss_bcphi(c,j+snl(c))
             mbc_pho(c,j) = mss_bcpho(c,j+snl(c))
             moc_phi(c,j) = mss_ocphi(c,j+snl(c))
             moc_pho(c,j) = mss_ocpho(c,j+snl(c))
             mdst1(c,j)   = mss_dst1(c,j+snl(c))
             mdst2(c,j)   = mss_dst2(c,j+snl(c))
             mdst3(c,j)   = mss_dst3(c,j+snl(c))
             mdst4(c,j)   = mss_dst4(c,j+snl(c))
             rds(c,j)     = snw_rds(c,j+snl(c))

          end if
       end do
    end do

    do fc = 1, num_snowc
       c = filter_snowc(fc)

       msno = abs(snl(c))

       if (msno == 1) then
          ! Specify a new snow layer
          if (dzsno(c,1) > 0.03_r8) then
             msno = 2
             dzsno(c,1) = dzsno(c,1)/2._r8
             swice(c,1) = swice(c,1)/2._r8
             swliq(c,1) = swliq(c,1)/2._r8
             dzsno(c,2) = dzsno(c,1)
             swice(c,2) = swice(c,1)
             swliq(c,2) = swliq(c,1)
             tsno(c,2)  = tsno(c,1)
             
             mbc_phi(c,1) = mbc_phi(c,1)/2._r8
             mbc_phi(c,2) = mbc_phi(c,1)
             mbc_pho(c,1) = mbc_pho(c,1)/2._r8
             mbc_pho(c,2) = mbc_pho(c,1)
             moc_phi(c,1) = moc_phi(c,1)/2._r8
             moc_phi(c,2) = moc_phi(c,1)
             moc_pho(c,1) = moc_pho(c,1)/2._r8
             moc_pho(c,2) = moc_pho(c,1)
             mdst1(c,1) = mdst1(c,1)/2._r8
             mdst1(c,2) = mdst1(c,1)
             mdst2(c,1) = mdst2(c,1)/2._r8
             mdst2(c,2) = mdst2(c,1)
             mdst3(c,1) = mdst3(c,1)/2._r8
             mdst3(c,2) = mdst3(c,1)
             mdst4(c,1) = mdst4(c,1)/2._r8
             mdst4(c,2) = mdst4(c,1)
             rds(c,2) = rds(c,1)

          end if
       end if

       if (msno > 1) then
          if (dzsno(c,1) > 0.02_r8) then
             drr = dzsno(c,1) - 0.02_r8
             propor = drr/dzsno(c,1)
             zwice = propor*swice(c,1)
             zwliq = propor*swliq(c,1)

             zmbc_phi = propor*mbc_phi(c,1)
             zmbc_pho = propor*mbc_pho(c,1)
             zmoc_phi = propor*moc_phi(c,1)
             zmoc_pho = propor*moc_pho(c,1)
             zmdst1 = propor*mdst1(c,1)
             zmdst2 = propor*mdst2(c,1)
             zmdst3 = propor*mdst3(c,1)
             zmdst4 = propor*mdst4(c,1)

             propor = 0.02_r8/dzsno(c,1)
             swice(c,1) = propor*swice(c,1)
             swliq(c,1) = propor*swliq(c,1)

             mbc_phi(c,1) = propor*mbc_phi(c,1)
             mbc_pho(c,1) = propor*mbc_pho(c,1)
             moc_phi(c,1) = propor*moc_phi(c,1)
             moc_pho(c,1) = propor*moc_pho(c,1)
             mdst1(c,1) = propor*mdst1(c,1)
             mdst2(c,1) = propor*mdst2(c,1)
             mdst3(c,1) = propor*mdst3(c,1)
             mdst4(c,1) = propor*mdst4(c,1)

             dzsno(c,1) = 0.02_r8

             mbc_phi(c,2) = mbc_phi(c,2)+zmbc_phi  ! (combo)
             mbc_pho(c,2) = mbc_pho(c,2)+zmbc_pho  ! (combo)
             moc_phi(c,2) = moc_phi(c,2)+zmoc_phi  ! (combo)
             moc_pho(c,2) = moc_pho(c,2)+zmoc_pho  ! (combo)
             mdst1(c,2) = mdst1(c,2)+zmdst1  ! (combo)
             mdst2(c,2) = mdst2(c,2)+zmdst2  ! (combo)
             mdst3(c,2) = mdst3(c,2)+zmdst3  ! (combo)
             mdst4(c,2) = mdst4(c,2)+zmdst4  ! (combo)
             rds(c,2) = rds(c,1) ! (combo)

             call Combo (dzsno(c,2), swliq(c,2), swice(c,2), tsno(c,2), drr, &
                  zwliq, zwice, tsno(c,1))

             ! Subdivide a new layer
             if (msno <= 2 .and. dzsno(c,2) > 0.07_r8) then
                msno = 3
                dtdz = (tsno(c,1) - tsno(c,2))/((dzsno(c,1)+dzsno(c,2))/2._r8) 
                dzsno(c,2) = dzsno(c,2)/2._r8
                swice(c,2) = swice(c,2)/2._r8
                swliq(c,2) = swliq(c,2)/2._r8
                dzsno(c,3) = dzsno(c,2)
                swice(c,3) = swice(c,2)
                swliq(c,3) = swliq(c,2)
                tsno(c,3) = tsno(c,2) - dtdz*dzsno(c,2)/2._r8
                if (tsno(c,3) >= tfrz) then 
                   tsno(c,3)  = tsno(c,2)
                else
                   tsno(c,2) = tsno(c,2) + dtdz*dzsno(c,2)/2._r8 
                endif

                mbc_phi(c,2) = mbc_phi(c,2)/2._r8
                mbc_phi(c,3) = mbc_phi(c,2)
                mbc_pho(c,2) = mbc_pho(c,2)/2._r8
                mbc_pho(c,3) = mbc_pho(c,2)
                moc_phi(c,2) = moc_phi(c,2)/2._r8
                moc_phi(c,3) = moc_phi(c,2)
                moc_pho(c,2) = moc_pho(c,2)/2._r8
                moc_pho(c,3) = moc_pho(c,2)
                mdst1(c,2) = mdst1(c,2)/2._r8
                mdst1(c,3) = mdst1(c,2)
                mdst2(c,2) = mdst2(c,2)/2._r8
                mdst2(c,3) = mdst2(c,2)
                mdst3(c,2) = mdst3(c,2)/2._r8
                mdst3(c,3) = mdst3(c,2)
                mdst4(c,2) = mdst4(c,2)/2._r8
                mdst4(c,3) = mdst4(c,2)
                rds(c,3) = rds(c,2)

             end if
          end if
       end if

       if (msno > 2) then
          if (dzsno(c,2) > 0.05_r8) then
             drr = dzsno(c,2) - 0.05_r8
             propor = drr/dzsno(c,2)
             zwice = propor*swice(c,2)
             zwliq = propor*swliq(c,2)
             
             zmbc_phi = propor*mbc_phi(c,2)
             zmbc_pho = propor*mbc_pho(c,2)
             zmoc_phi = propor*moc_phi(c,2)
             zmoc_pho = propor*moc_pho(c,2)
             zmdst1 = propor*mdst1(c,2)
             zmdst2 = propor*mdst2(c,2)
             zmdst3 = propor*mdst3(c,2)
             zmdst4 = propor*mdst4(c,2)

             propor = 0.05_r8/dzsno(c,2)
             swice(c,2) = propor*swice(c,2)
             swliq(c,2) = propor*swliq(c,2)

             mbc_phi(c,2) = propor*mbc_phi(c,2)
             mbc_pho(c,2) = propor*mbc_pho(c,2)
             moc_phi(c,2) = propor*moc_phi(c,2)
             moc_pho(c,2) = propor*moc_pho(c,2)
             mdst1(c,2) = propor*mdst1(c,2)
             mdst2(c,2) = propor*mdst2(c,2)
             mdst3(c,2) = propor*mdst3(c,2)
             mdst4(c,2) = propor*mdst4(c,2)

             dzsno(c,2) = 0.05_r8

             mbc_phi(c,3) = mbc_phi(c,3)+zmbc_phi  ! (combo)
             mbc_pho(c,3) = mbc_pho(c,3)+zmbc_pho  ! (combo)
             moc_phi(c,3) = moc_phi(c,3)+zmoc_phi  ! (combo)
             moc_pho(c,3) = moc_pho(c,3)+zmoc_pho  ! (combo)
             mdst1(c,3) = mdst1(c,3)+zmdst1  ! (combo)
             mdst2(c,3) = mdst2(c,3)+zmdst2  ! (combo)
             mdst3(c,3) = mdst3(c,3)+zmdst3  ! (combo)
             mdst4(c,3) = mdst4(c,3)+zmdst4  ! (combo)
             rds(c,3) = rds(c,2) ! (combo)

             call Combo (dzsno(c,3), swliq(c,3), swice(c,3), tsno(c,3), drr, &
                  zwliq, zwice, tsno(c,2))

             ! Subdivided a new layer
             if (msno <= 3 .and. dzsno(c,3) > 0.18_r8) then
                msno =  4
                dtdz = (tsno(c,2) - tsno(c,3))/((dzsno(c,2)+dzsno(c,3))/2._r8) 
                dzsno(c,3) = dzsno(c,3)/2._r8
                swice(c,3) = swice(c,3)/2._r8
                swliq(c,3) = swliq(c,3)/2._r8
                dzsno(c,4) = dzsno(c,3)
                swice(c,4) = swice(c,3)
                swliq(c,4) = swliq(c,3)
                tsno(c,4) = tsno(c,3) - dtdz*dzsno(c,3)/2._r8
                if (tsno(c,4) >= tfrz) then 
                   tsno(c,4)  = tsno(c,3)
                else
                   tsno(c,3) = tsno(c,3) + dtdz*dzsno(c,3)/2._r8 
                endif
                
                mbc_phi(c,3) = mbc_phi(c,3)/2._r8
                mbc_phi(c,4) = mbc_phi(c,3)
                mbc_pho(c,3) = mbc_pho(c,3)/2._r8
                mbc_pho(c,4) = mbc_pho(c,3)
                moc_phi(c,3) = moc_phi(c,3)/2._r8
                moc_phi(c,4) = moc_phi(c,3)
                moc_pho(c,3) = moc_pho(c,3)/2._r8
                moc_pho(c,4) = moc_pho(c,3)
                mdst1(c,3) = mdst1(c,3)/2._r8
                mdst1(c,4) = mdst1(c,3)
                mdst2(c,3) = mdst2(c,3)/2._r8
                mdst2(c,4) = mdst2(c,3)
                mdst3(c,3) = mdst3(c,3)/2._r8
                mdst3(c,4) = mdst3(c,3)
                mdst4(c,3) = mdst4(c,3)/2._r8
                mdst4(c,4) = mdst4(c,3)
                rds(c,4) = rds(c,3)

             end if
          end if
       end if

       if (msno > 3) then
          if (dzsno(c,3) > 0.11_r8) then
             drr = dzsno(c,3) - 0.11_r8
             propor = drr/dzsno(c,3)
             zwice = propor*swice(c,3)
             zwliq = propor*swliq(c,3)
             
             zmbc_phi = propor*mbc_phi(c,3)
             zmbc_pho = propor*mbc_pho(c,3)
             zmoc_phi = propor*moc_phi(c,3)
             zmoc_pho = propor*moc_pho(c,3)
             zmdst1 = propor*mdst1(c,3)
             zmdst2 = propor*mdst2(c,3)
             zmdst3 = propor*mdst3(c,3)
             zmdst4 = propor*mdst4(c,3)

             propor = 0.11_r8/dzsno(c,3)
             swice(c,3) = propor*swice(c,3)
             swliq(c,3) = propor*swliq(c,3)

             mbc_phi(c,3) = propor*mbc_phi(c,3)
             mbc_pho(c,3) = propor*mbc_pho(c,3)
             moc_phi(c,3) = propor*moc_phi(c,3)
             moc_pho(c,3) = propor*moc_pho(c,3)
             mdst1(c,3) = propor*mdst1(c,3)
             mdst2(c,3) = propor*mdst2(c,3)
             mdst3(c,3) = propor*mdst3(c,3)
             mdst4(c,3) = propor*mdst4(c,3)

             dzsno(c,3) = 0.11_r8

             mbc_phi(c,4) = mbc_phi(c,4)+zmbc_phi  ! (combo)
             mbc_pho(c,4) = mbc_pho(c,4)+zmbc_pho  ! (combo)
             moc_phi(c,4) = moc_phi(c,4)+zmoc_phi  ! (combo)
             moc_pho(c,4) = moc_pho(c,4)+zmoc_pho  ! (combo)
             mdst1(c,4) = mdst1(c,4)+zmdst1  ! (combo)
             mdst2(c,4) = mdst2(c,4)+zmdst2  ! (combo)
             mdst3(c,4) = mdst3(c,4)+zmdst3  ! (combo)
             mdst4(c,4) = mdst4(c,4)+zmdst4  ! (combo)
             rds(c,4) = rds(c,3) ! (combo)

             call Combo (dzsno(c,4), swliq(c,4), swice(c,4), tsno(c,4), drr, &
                  zwliq, zwice, tsno(c,3))

             ! Subdivided a new layer
             if (msno <= 4 .and. dzsno(c,4) > 0.41_r8) then
                msno = 5
                dtdz = (tsno(c,3) - tsno(c,4))/((dzsno(c,3)+dzsno(c,4))/2._r8) 
                dzsno(c,4) = dzsno(c,4)/2._r8
                swice(c,4) = swice(c,4)/2._r8
                swliq(c,4) = swliq(c,4)/2._r8
                dzsno(c,5) = dzsno(c,4)
                swice(c,5) = swice(c,4)
                swliq(c,5) = swliq(c,4)
                tsno(c,5) = tsno(c,4) - dtdz*dzsno(c,4)/2._r8 
                if (tsno(c,5) >= tfrz) then 
                   tsno(c,5)  = tsno(c,4)
                else
                   tsno(c,4) = tsno(c,4) + dtdz*dzsno(c,4)/2._r8 
                endif

                mbc_phi(c,4) = mbc_phi(c,4)/2._r8
                mbc_phi(c,5) = mbc_phi(c,4)
                mbc_pho(c,4) = mbc_pho(c,4)/2._r8
                mbc_pho(c,5) = mbc_pho(c,4)              
                moc_phi(c,4) = moc_phi(c,4)/2._r8
                moc_phi(c,5) = moc_phi(c,4)
                moc_pho(c,4) = moc_pho(c,4)/2._r8
                moc_pho(c,5) = moc_pho(c,4)
                mdst1(c,4) = mdst1(c,4)/2._r8
                mdst1(c,5) = mdst1(c,4)
                mdst2(c,4) = mdst2(c,4)/2._r8
                mdst2(c,5) = mdst2(c,4)
                mdst3(c,4) = mdst3(c,4)/2._r8
                mdst3(c,5) = mdst3(c,4)
                mdst4(c,4) = mdst4(c,4)/2._r8
                mdst4(c,5) = mdst4(c,4)
                rds(c,5) = rds(c,4)

             end if
          end if
       end if

       if (msno > 4) then
          if (dzsno(c,4) > 0.23_r8) then
             drr = dzsno(c,4) - 0.23_r8
             propor = drr/dzsno(c,4)
             zwice = propor*swice(c,4)
             zwliq = propor*swliq(c,4)
             
             zmbc_phi = propor*mbc_phi(c,4)
             zmbc_pho = propor*mbc_pho(c,4)
             zmoc_phi = propor*moc_phi(c,4)
             zmoc_pho = propor*moc_pho(c,4)
             zmdst1 = propor*mdst1(c,4)
             zmdst2 = propor*mdst2(c,4)
             zmdst3 = propor*mdst3(c,4)
             zmdst4 = propor*mdst4(c,4)

             propor = 0.23_r8/dzsno(c,4)
             swice(c,4) = propor*swice(c,4)
             swliq(c,4) = propor*swliq(c,4)

             mbc_phi(c,4) = propor*mbc_phi(c,4)
             mbc_pho(c,4) = propor*mbc_pho(c,4)
             moc_phi(c,4) = propor*moc_phi(c,4)
             moc_pho(c,4) = propor*moc_pho(c,4)
             mdst1(c,4) = propor*mdst1(c,4)
             mdst2(c,4) = propor*mdst2(c,4)
             mdst3(c,4) = propor*mdst3(c,4)
             mdst4(c,4) = propor*mdst4(c,4)

             dzsno(c,4) = 0.23_r8

             mbc_phi(c,5) = mbc_phi(c,5)+zmbc_phi  ! (combo)
             mbc_pho(c,5) = mbc_pho(c,5)+zmbc_pho  ! (combo)
             moc_phi(c,5) = moc_phi(c,5)+zmoc_phi  ! (combo)
             moc_pho(c,5) = moc_pho(c,5)+zmoc_pho  ! (combo)
             mdst1(c,5) = mdst1(c,5)+zmdst1  ! (combo)
             mdst2(c,5) = mdst2(c,5)+zmdst2  ! (combo)
             mdst3(c,5) = mdst3(c,5)+zmdst3  ! (combo)
             mdst4(c,5) = mdst4(c,5)+zmdst4  ! (combo)
             rds(c,5) = rds(c,4) ! (combo)

             call Combo (dzsno(c,5), swliq(c,5), swice(c,5), tsno(c,5), drr, &
                  zwliq, zwice, tsno(c,4))
          end if
       end if

       snl(c) = -msno

    end do

    do j = -nlevsno+1,0
       do fc = 1, num_snowc
          c = filter_snowc(fc)
          if (j >= snl(c)+1) then
             dz(c,j)         = dzsno(c,j-snl(c))
             h2osoi_ice(c,j) = swice(c,j-snl(c))
             h2osoi_liq(c,j) = swliq(c,j-snl(c))
             t_soisno(c,j)   = tsno(c,j-snl(c))

             mss_bcphi(c,j)   = mbc_phi(c,j-snl(c))
             mss_bcpho(c,j)   = mbc_pho(c,j-snl(c))
             mss_ocphi(c,j)   = moc_phi(c,j-snl(c))
             mss_ocpho(c,j)   = moc_pho(c,j-snl(c))
             mss_dst1(c,j)    = mdst1(c,j-snl(c))
             mss_dst2(c,j)    = mdst2(c,j-snl(c))
             mss_dst3(c,j)    = mdst3(c,j-snl(c))
             mss_dst4(c,j)    = mdst4(c,j-snl(c))
             snw_rds(c,j)     = rds(c,j-snl(c))

          end if
       end do
    end do

    do j = 0, -nlevsno+1, -1
       do fc = 1, num_snowc
          c = filter_snowc(fc)
          if (j >= snl(c)+1) then
             z(c,j)    = zi(c,j) - 0.5_r8*dz(c,j)
             zi(c,j-1) = zi(c,j) - dz(c,j)
          end if
       end do
    end do

  end subroutine DivideSnowLayers

!-----------------------------------------------------------------------
!BOP
!
! !IROUTINE: Combo
!
! !INTERFACE:
  subroutine Combo(dz,  wliq,  wice, t, dz2, wliq2, wice2, t2)
!
! !DESCRIPTION:
! Combines two elements and returns the following combined
! variables: dz, t, wliq, wice.
! The combined temperature is based on the equation:
! the sum of the enthalpies of the two elements =
! that of the combined element.
!
! !USES:
    use clm_varcon,  only : cpice, cpliq, tfrz, hfus
!
! !ARGUMENTS:
    implicit none
    real(r8), intent(in)    :: dz2   ! nodal thickness of 2 elements being combined [m]
    real(r8), intent(in)    :: wliq2 ! liquid water of element 2 [kg/m2]
    real(r8), intent(in)    :: wice2 ! ice of element 2 [kg/m2]
    real(r8), intent(in)    :: t2    ! nodal temperature of element 2 [K]
    real(r8), intent(inout) :: dz    ! nodal thickness of 1 elements being combined [m]
    real(r8), intent(inout) :: wliq  ! liquid water of element 1
    real(r8), intent(inout) :: wice  ! ice of element 1 [kg/m2]
    real(r8), intent(inout) :: t     ! nodel temperature of elment 1 [K]
!
! !CALLED FROM:
! subroutine CombineSnowLayers in this module
! subroutine DivideSnowLayers in this module
!
! !REVISION HISTORY:
! 15 September 1999: Yongjiu Dai; Initial code
! 15 December 1999:  Paul Houser and Jon Radakovich; F90 Revision
!
!
! !LOCAL VARIABLES:
!EOP
!
    real(r8) :: dzc   ! Total thickness of nodes 1 and 2 (dzc=dz+dz2).
    real(r8) :: wliqc ! Combined liquid water [kg/m2]
    real(r8) :: wicec ! Combined ice [kg/m2]
    real(r8) :: tc    ! Combined node temperature [K]
    real(r8) :: h     ! enthalpy of element 1 [J/m2]
    real(r8) :: h2    ! enthalpy of element 2 [J/m2]
    real(r8) :: hc    ! temporary
!-----------------------------------------------------------------------

    dzc = dz+dz2
    wicec = (wice+wice2)
    wliqc = (wliq+wliq2)
    h = (cpice*wice+cpliq*wliq) * (t-tfrz)+hfus*wliq
    h2= (cpice*wice2+cpliq*wliq2) * (t2-tfrz)+hfus*wliq2

    hc = h + h2
    tc = tfrz + (hc - hfus*wliqc) / (cpice*wicec + cpliq*wliqc)

    dz = dzc
    wice = wicec
    wliq = wliqc
    t = tc

  end subroutine Combo

!-----------------------------------------------------------------------
!BOP
!
! !IROUTINE: BuildSnowFilter
!
! !INTERFACE:
  subroutine BuildSnowFilter(lbc, ubc, num_nolakec, filter_nolakec, &
                             num_snowc, filter_snowc, &
                             num_nosnowc, filter_nosnowc)
!
! !DESCRIPTION:
! Constructs snow filter for use in vectorized loops for snow hydrology.
!
! !USES:
    use clmtype
!
! !ARGUMENTS:
    implicit none
    integer, intent(in)  :: lbc, ubc                    ! column bounds
    integer, intent(in)  :: num_nolakec                 ! number of column non-lake points in column filter
    integer, intent(in)  :: filter_nolakec(ubc-lbc+1)   ! column filter for non-lake points
    integer, intent(out) :: num_snowc                   ! number of column snow points in column filter
    integer, intent(out) :: filter_snowc(ubc-lbc+1)     ! column filter for snow points
    integer, intent(out) :: num_nosnowc                 ! number of column non-snow points in column filter
    integer, intent(out) :: filter_nosnowc(ubc-lbc+1)   ! column filter for non-snow points
!
! !CALLED FROM:
! subroutine Hydrology2 in Hydrology2Mod
! subroutine CombineSnowLayers in this module
!
! !REVISION HISTORY:
! 2003 July 31: Forrest Hoffman
!
! !LOCAL VARIABLES:
! local pointers to implicit in arguments
    integer , pointer :: snl(:)                        ! number of snow layers
!
!
! !OTHER LOCAL VARIABLES:
!EOP
    integer  :: fc, c
!-----------------------------------------------------------------------

    ! Assign local pointers to derived subtype components (column-level)

    snl => clm3%g%l%c%cps%snl

    ! Build snow/no-snow filters for other subroutines

    num_snowc = 0
    num_nosnowc = 0
    do fc = 1, num_nolakec
       c = filter_nolakec(fc)
       if (snl(c) < 0) then
          num_snowc = num_snowc + 1
          filter_snowc(num_snowc) = c
       else
          num_nosnowc = num_nosnowc + 1
          filter_nosnowc(num_nosnowc) = c
       end if
    end do

  end subroutine BuildSnowFilter

end module SnowHydrologyMod
